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ABSTRACT 

We study the torque on low mass protoplanets on fixed circular orbits, embedded in 
a protoplanetary disc in the isothermal limit. We consider a wide range of surface 
density distributions including cases where the surface density increases smoothly 
outwards. We perform both linear disc response calculations and non linear numerical 
simulations. We consider a large range of viscosities, including the inviscid limit, as 
well as a range of protoplanet mass ratios, with special emphasis on the coorbital 
region and the corotation torque acting between disc and protoplanet. 

For low mass protoplanets and large viscosity the corotation torque behaves as 
expected from linear theory. However, when the viscosity becomes small enough to 
enable horseshoe turns to occur, the linear corotation torque exists only temporarily 
after insertion of a planet into the disc, being replaced by the horseshoe drag first 
discussed by Ward. This happens after a time that is equal to the horseshoe libra- 
tion period reduced by a factor amounting to about twice the disc aspect ratio. This 
torque scales with the radial gradient of specific vorticity, as does the linear torque, 
but we find it to be many times larger. If the viscosity is large enough for viscous 
diffusion across the coorbital region to occur within a libration period, we find that 
the horseshoe drag may be sustained. If not, the corotation torque saturates leaving 
only the linear Lindblad torques. As the magnitude of the non linear coorbital torque 
(horseshoe drag) is always found to be larger than the linear torque, we find that 
the sign of the total torque may change even for for mildly positive surface density 
gradients. In combination with a kinematic viscosity large enough to keep the torque 
from saturating, strong sustained deviations from linear theory and outward or stalled 
migration may occur in such cases. 

Key words: planetary systems: formation - planets and satellites: formation. 



1 INTRODUCTION 

A planet embedded in a gaseous disc is subject to a net 
torque that gives rise to orb ital evolution. Since th e discov- 
ery of the first hot Jupiter (jMavor fe Queloall995l ). planet 
migration due to interaction with the protoplanetary disc 
has become a necessary ingredient of planet formation the- 
ory. A considerable amount of analytical and numerical work 
has gone into understanding t he disc-planet interact ions 
leading to planet migration (see iPapaloizou et al.ll2007l . for 
an overview). 

Three modes of migration can be distinguished, in 
most cases leading to migration towards the central star. 
High mass planets, for which the Hill sphere is larger 
than the disc scale height, open up deep gaps in the 
disc, after which migration proceeds on a viscous time 
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scale (|Lin fc Papaloi zou 1986). For standard disc param- 
eters, planets mo re massive than Jupiter migrate in this 
Type II regime i|Crida fc Morbidellil 120071 ). Intermediate 
mass planets, comparable to Saturn, embedded in mas- 
sive discs may undergo ru naway or Type III migration 
IjMasset fc Papaloizoull2003|h which may be directed inward 
as well as outward l|Pepliriski et alj|200cf ). In this paper, we 
focus on low mass planets, with masses typical up to a few 
times the mass of the Earth (Mq). Low mass planets excite 
linear wav es in the dis c, the action of which leads to Type I 
migration (Ward 1993). Linear, semi-analytic al calculations 
have resulted in a widely used torque formula (|Tanaka et al.l 
l2002i . hereafter TTW02) for isothermal discs. 

The total torque on an embedded planet can be decom- 
posed in torques due to waves, which are generated at Lind- 
blad resonances and propagate away from the planet, and 
corotation torques, generated near the orbit of the planet, 
where material on average corotates with the planet (see 
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iGoldreich fc Tremainelll979l ). For low mass planets, the wave 
torque is usually thought to dominate the total torque, and, 
being relatively insensitive to background gradients, lead to 
inward migration. 

Although Type I migration was thought to be math- 
ematically well-understood, it nevertheless posed a serious 
problem for planet form ation. Applying the torque fo rmula 
from TTW02 (see also iKorvcanskv fc Pollack! 1 1991 here- 
after KP93) to planets of a few M® embedded in a typi- 
cal disc resulted in inward migration time scales that are 
much shorter than the lifetime of the disc. In other words, 
these planets would all migr ate to very small orbital radii, 
or even into the central star (Ward 1997). Since gaseous gi- 
ant planets are thought to form around a solid core that 
is in this mass range, Type I migration theory essentially 
predicts that there should be no planets at large radii. This 
has led to investigations on how to slow down or stop Type I 
migration, for example through the action o f magnetic fields 
(|Terauemll2003l ; lNelson fc Papaloizoull2004 ) or sharp surface 
density gradients ( Masset et al.l 1200*61 ). 

More recently, efforts have been made to relax the 
isothermal assumption that was made in previous works, 
and to properly account for the energy budget of the disc. 
Since radiation is the dominant cooling agent, radiation- 
hydrodynamical simulations are needed. The outcome of 
these si mulations were surprising: iPaardekooper fc Mellerrial 
(2006a) found that the torque on the planet was positive, 
leading to outward migration, whenever the opacity of the 
disc at the location of the planet was high enoug h to m ake 
cooling inefficient. In IPaardekooper fc Mellemal (|2008l ) it 
was recognised that this positive torque was a corotation ef- 
fect, related to a radial en t ropy gradient in the unperturbed 
disc. iBaruteau fc Massed (|2008l ). through a linear analysis 
of the corotation torque, suggested that the linear corota- 
tion torque in the presence of an entropy gradient can be 
stron g enough to overcome the ne g ative wave torque. How- 
ever, IPaardeko oper fc Papaloizou! (|2008l ) showed that this 
linear contribution is small, and that a genuinely non-linear 
effect is responsible for the change of sign of the total torque. 

In this paper, we take a step back and reanalyse the 
isothermal case, in a two-dimensional set-up. The simplic- 
ity of this model allows us to perform numerical simulations 
with high enough resolution investigate in detail non-linear 
effects on the corotation torque and to run them for long 
enough to study its possible saturation. We will show that 
linear theory only applies for short times on the order of 
a few orbits when a protoplanet is inserted into a disc for 
which the viscosity is not too large and there is a non-zero 
corotation torque. We also show t hat it is a non -linear effect 
associated with horseshoe bends (|Wardlll99ll ) that results 
in a departure from linear theory. This departure is found 
to result in torques that can have a much larger magnitude 
than the linear ones. As the torque scales with the gradient 
of specific vorticity, this may produce noticeable effects that 
may even lead to stalled or outward migration, when this 
gradient is relatively large as for example occurs when the 
surface density increases gently outwards. For these effects 
to be sustained the viscosity must be large enough to pre- 
vent torque saturation. Viscosities comparable to those often 
assumed in protoplanetary disc modelling are found to be 
large enough to enable these non linear corotation torques 
to be sustained for protoplanets in the Earth mass range. 



The plan of the paper is as follows. We start in Sec- 
tion [2] by reviewing the model used throughout the pa- 
per. In Section [3] we review and discuss linear corotation 
torque estimates as well as torque estimates based on the 
horseshoe drag experienced by material executing horseshoe 
turns. These are shown to be distinct phenomena with dif- 
ferent dependencies on the physical variables, the horseshoe 
drag being essentially non linear even though the associated 
torque scales in the same way as the linear one. We go on to 
perform detailed linear calculations in Section [4] and then 
compare the linear and non-linear (horseshoe drag) torques 
in Section [S] In Section HJ] we present the results of numerical 
hydrodynamical simulations, torques are also compared to 
linear and horseshoe drag estimates confirming the view out- 
lined above. The horseshoe drag is found to be significantly 
larger than the linear corotation torque and potentially a 
very important contributor to the total torque for disc sur- 
face densities that increase even mildly outwards. We go on 
to consider the long term torque evolution and saturation as 
a function of viscosity finding corotation torque saturation 
in low viscosity cases and sustained corotation torques when 
the viscosity is large enough to resupply angular momentum 
to the coorbital region. We discuss our results further in Sec- 
tion [7] Finally we present a short summary, together with 
some concluding remarks, in Section [8] In addition we give 
an analysis and discussion of the time development of the 
linear corotation torque acting on a planet after immersion 
into a disc, demonstrating that the characteristic time is the 
orbital period as was found in the numerical simulations. 



2 BASIC EQUATIONS 

The evolution of a gaseous disc is governed by the Navier- 
Stokes equations. We will work in a cylindrical coordinate 
polar frame (r,ip), centred on the central star, and we in- 
tegrate the equations of motion vertically to obtain a two- 
dimensional problem. We are then left with the continuity 
equation and two equations of motion: 



dv 

dt 



+ (v • V)v = 



dt 

Vp 
~E~ 



+ V • Ev = 



V$ + f vis 



(1) 

(2) 



where £ is the surface density, v = («, rQ) T is the velocity, 
p is the vertically integrated pressure, $ is the gravitational 
potential and f v i sc represents the viscous force. We use a lo- 
cally isothermal equation of state, p = c^E, where the sound 
speed c s may be a prescribed function of radius. We take dj 
to be a power law with index —j3, which makes — f3 basically 
the power law index of the radial temperature profile. Writ- 
ing c s = HQk, where H is the vertical pressure scale height 
and J1k is the Keplerian angular velocity, we usually adopt 
either a sound speed that gives rise to a constant aspect 
ratio h = H/r (or, equivalently, (3 = 1), or a purely isother- 
mal disc with constant c s (/3 = 0). In the latter case, we 
usually quote the aspect ratio at the location of the planet, 
hp = c s /(r p f2 p ). The initial, or for linear calculations the 
background, surface density is taken to be a power law with 
index —a. Thus £ = E p (r/r p ) _a , where E p is the surface 
density at the location of the planet. The gravitational po- 
tential $ is the sum of the potential due to the central star, 
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the planet's potential $ p , and an indirect term that arises 
due the the acceleration of the coordinate frame centred on 
the star. For $ p , we use a softened point mass potential: 



GM P 



, (3) 

^Jr 2 — 2rr p cos((p — (p p ) + r p + b 2 r p 

where G is the gravitational constant, M p is the mass of the 
planet, and (r p ,ip p ) are the coordinates of the planet. We 
will also use q to indicate the mass ratio M p /M«, where M* 
is the mass of the central star. The softening parameter b 
should be a sizeable fraction of h, in order to approximate 
the result of appropriate vertical averaging of the potential. 

The ex act form of the viscou s terms can be found else- 
where (e.g. iD'Aneelo et al.1 f2002\ We take the kinematic 
viscosity v to be a power law in radius, such that the initial 
surface density profile is a stationary solution. It is easy to 
see that the required viscosity law is v cx r a_1,/2 . This way, 
we can neglect a global radial velocity field in the model, 
which greatly simplifies the analysis. Values quoted in the 
text refer to v at the location of the planet, and will be in 
units of r p Q p , where f2 p is the angular velocity of the planet. 



3 COROTATION TORQUES 

The wave (or Lindblad) torque exerted by a low-mass 
planet on a gaseous disc i s relatively well-understood 
(|Goldreich fe Tremaind Il979l , TTW02). When the Hill 
sphere of the planet is much smaller than the scale height 
of the disc, the density waves are excited at Lindblad res- 
onances are linear and the resulting torque can be calcu- 
lated by performing a linear response calculation and then 
summing the contributions arising from individual Fourier 
components (TTW02). This torque usually leads to inward 
migration, and has been thought to dominate in the Type I 
regime of low mass planets. There are two approaches that 
can be followed in order to obtain corotation torques. The 
first, which we will refer to as the linear estimate, is based 
on a linear response calculation. The second is a fundamen- 
tally different approach based on a direct torque calculation 
made after making some assumptions about the trajectories 
of the disc fluid elements, which we refer to as a horseshoe 
drag calculation. The first approach applies at early times 
after a protoplanet is inserted into a disc as, at least for a 
low mass protoplanet, the response must first of all be linear. 
The second approach applies at later times after the stream- 
line pattern has adjusted to the presence of the planet. 



3.1 Formula for the linear corotation torque for 
each azimuthal mode number 

In order to obtain linear estimates for corotation torques, 
the equations are linearised (see Section [4} and Fourier- 
decomposed in azimuth. The resulting sing le second-order 
equation (see iGoldreich fc Tremainel [l979) governing the 
disc response can be solved to find the corotation torque 
acting on the planet. Using an approximation scheme that 
assumes the perturbing potential v aries on a scale signifi- 
cantly longer then the scale height, IGoldreich fc Tremaind 
( 1979) find this torque given by 



where <&' m is the m-th Fourier component of the planet's 
potential, m is the azimuthal mode number and ui is the 
flow vorticity, being equal to Q/2 in a Keplerian disc. All 
quantities in equation Q should be evaluated at corota- 
tion (r = r p ), which makes the total torque proportional 
to the radial gradient in specific vorticity, or vortensity, in 
the unperturbed flow there. Of course unless the softening 
parameter is rather large, the perturbing potential does not 
vary slowly on the corotation circle and so @ cannot be 
used as described above. However, it becomes valid if the 
perturbing potential <&' m is replaced by the generalised po- 
tential obtained by adding the enthalpy perturbation to it 
(see e.g. lLai fc Zhandl2006l . and the appendix). But then the 
linear response equations need to be solved in order to deter- 
mine the enthalpy perturbation in order to evaluate torques 
using the modified form of Q, which are nonetheless still 
proportional to the gradient of specific vorticity. 

3.2 The total linear corotation torque in the limit 
of zero softening 

To obtain the total torque, on e needs to s um the contribu- 
tions from all values of m (see IWardlll989l ). For a Keplerian 
disc with zero softening, the total torque has been calcu- 
lated, by finding the linear response numerically, to be given 
by (see TTW02) 

„2 



1.36 



(5) 



r — 

J- cm — 



TO7T d 

2dQ./dr dr \u 



(4) 



For a 3D disc it is found that the same expression holds but 
with the numerical coefficient being replaced by 0.632. 

An important issue is the time required to develop the 
linear corotation torque. Being linear this should not de- 
pend on the mass of the protoplanet but only on intrinsic 
disc parameters. Furthermore when a protoplanet is inserted 
into a disc, there has to be an initial linear phase and for 
sufficiently low protoplanet mass, the full linear response 
should develop. As it is somewhat involved, we relegate the 
discussion of the time development of the linear corotation 
torque to the Appendix. From this discussion, the charac- 
teristic development time is expected to be on the order of 
the orbital period. We now go on to consider the subsequent 
development of the corotational flow. 

3.3 The horseshoe drag 

IWardl l|l99lf ) derived a formula for the corotation torque us- 
ing a fundamentally different approach. The argument was 
based on the expected form of the gas streamlines. These 
can be classified into four groups as viewed in a reference 
frame corotating with the planet (see Fig. [TJ. The first 
pass the protoplanet interior to the coorbital region and 
the second pass it exterior to the coorbital region. These 
groups extend to large distances from the protoplanet. The 
other two groups consist of material in the corotation re- 
gion on horseshoe orbits that execute turns close to the 
planet. The third group approach and leave the protoplanet 
from its leading side while the fourth do so from its trailing 
side. These groups are separated by two separatrices (see 
e.g.. iMasset et afll2006l ; iPaardekooper fc Papaloizoul [20081 . 
for more discussion of this aspect). For low mass protoplan- 
ets, most of the corotation torque is produced in a region 
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Figure 1. Schematic overview of the horseshoe region, with the 
planet denoted by the black circle. Solid grey lines indicate ap- 
proximate streamlines, with the separatrices in black. The sepa- 
ratices divide the region into four parts, labelled by Roman nu- 
merals: I and II, the inner and the outer disc, respectively, and 
III and IV the leading and trailing part of the horseshoe region, 
respectively. Two horizontal dashed lines indicate the boundaries 
of the torque calculation, Ti and I^. Note that in this figure, the 
location of Ti and T2 is chosen arbitrarily. In the model, they are 
chosen to be far enough from the planet so that the corotation 
torque is determined within. 



close to the planet with a length scale expected to be the 
larger of br p or H. It is important to note that it takes 
some time for the streamline structure described above to 
develop on this scale after insertion of a protoplanet into 
a disc. As this structure represents a finite deviation from 
the initial form, this time depends on the mass of the proto- 
planet. The calculation of the horseshoe drag applies to the 
situation when this streamline structure has developed on 
this scale. It is important to note that the time required is 
shorter than that required to develop it on a scale compara- 
ble to r p which is when issues of torque saturation need to 

be considered. 

Following Ward] (|l99ll ) we consider the torque produced 
by material on streamlines undergoing horseshoe turns. We 
consider a region 1Z interior to the two separatrices sepa- 
rating these from the first and fourth group of streamlines 
that is bounded by two lines of constant ip, Fi and T2 on 
the trailing and leading sides of the protoplanet respectively 
(see Fig.[T]). These boundaries are supposed to be sufficiently 
far from the protoplanet that the corotation torque is deter- 
mined within. Assuming a steady state this torque may be 
obtained by considering the conservation of angular momen- 
tum within 1Z written in the form 



Tc.hs — 



E(j -j p )(n-fl p )rdr 



(6) 



Here j = rv v is the specific angular momentum and j p is j 
evaluated at the or bital radius o f the protoplanet. 

We now follow IWardl (|l99lh and assume the streamline 
pattern is symmetric on the leading and trailing sides such 
that a streamline entering on T2 can be identified with a 
corresponding streamline leaving on Ti at the same radial 
location. But note that on account of the horseshoe turns, 
these streamlines will enter 7Z at different radii. Since, in 
the barotropic case, potential vorticity or vortensity E/lj 
is conserved along streamlines, and these streamlines are 
disconnected, this will differ on them. Accordingly we write 
the torque as 



Here 



L 



^ A y- ) u(j - j p )(Cl - Q p )rdr. 
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(7) 



(8) 



is the vortensity difference on the corresponding streamlines 
which have been assumed to enter 7Z at the same radial 



distance from the planet on opposite sides, 



Q p /2 is the 



vorticity at the planets orbital radius and we have assumed 
the vortensity to be an even function of radial distance from 
the protoplanet. 

The above integral is easily done if one adopts a first 
order Taylor expansion for the quantities in brackets and 



less w idth of the horseshoe region is x s . One obtains (|Wardl 
1 19911 ): 



(9) 



Note that as is apparent from the above discussion and 
that given in the Appendix and also the results of nu- 
merical simulations to be presented later, the horse- 
shoe drag, given by equation ©, occurs as a non 
linear e ffect that has no counterpart in linear theory 
(see also Paardekooper fc Papaloizoul [2008. in addition to 
iPaardekooper fc Papaloizoul 120091 ). We also note that nu- 
merical hydrodynamical calculations necessarily have b > 0, 
and, two-dimensional simulations, generally adopt b ~ h 
to account for three-dimensional effects in an approximate 
way. We will see that this strongly affects the width of the 
horseshoe region x s . A detailed analysis of the horsesho e 
region is presented in IPaardekooper fc Papaloizoul {2009 ). 
Here, we adopt a simple estimate for x s that has proved 
to be reasonable for smoothing le ngths comparable to h p 
l| Paardekooper fc Pap aloizou 2008): 



(10) 



In general, this dependence is to be expected when the sep- 
aratrix streamline passes through the location of the planet 
wit h b possibly being r eplaced by h p for small softening (see 
also lMasset et af]|2006l ). We will compare the horseshoe drag 
resulting from x s to the linear corotation torque in Section 
\5\ Note that, since r c ,hs is proportional to a;*, the horse- 
shoe drag would then be proportional to q 2 /h p in the small 
softening case and q 2 /b 2 in the large softening case , just as 
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r c ,n n . This means that the dependence of the total torque 
on q, b and h p would not be a good indication of linearity, 
since both the linear and the non-linear corotation torque 
scale in the same way. The only ways to distinguish these 
are through their magnitudes, and by their time evolution. 

The linear corotation torque is set up in approximately 
an orbital time scale (see the Appendix), similar to the wave 
torque. This means that, when the background state of the 
disc is stationary, any evolution in the torque after a few 
dynamical time scales is due to non-linear effects. The finite 
width of the horseshoe region gives rise to a libration time 
scale 

•nib = -j^-^p 1 ' ( n ) 

which is basically the time it takes for a fluid element at 
orbital radius r p (l + x s ) to complete two orbits in a frame 
corotating with the planet. Therefore, this is the time scale 
on which t he corotatio n torque will saturate due to phase 
mixing rsee lWard|[2007h . unl ess some form of viscosity oper- 
ates on smaller time scales (jMassct 2002). Note tha t satu- 
ration is a non-linear process (|Qgilvie fc Lubowl l2003) , since 
it hinges on the finite width of the corotation region. 

We will see that there are basically three time scales 
in the problem, two of which are closely related to non- 
linearity. First, there is the orbital time scale, on which all 
linear torques are set up. This is the shortest time scale in 
the problem. The longest time scale is rub, on which satu- 
ration operates. A third time scale governs the development 
of the horseshoe drag, w hich we will see takes a fraction of a 
libration time (see also iPaardekooper fc Papaloizou| [2008). 
and lies in between the time scale for the linear torque de- 
velopment and the saturation time scale. 



4 LINEAR CALCULATIONS 

Since we are interested in departures from linear theory, it 
is necessary to first firmly establish what linear theory pre- 
dicts. The 2D linear calculations of KP93 and TTW02 use a 
vanishingly small value for b, and these results are therefore 
not directly comparable to our hydrodynamical simulations. 
Our customised linear calculations, as briefly outlined be- 
low, with b comparable to h, can be directly compared to 
hydrodynamical simulations. 



4.1 Governing equations 

Linearising equations {TJ and with f v i sc = 0, and assum- 
ing a Fourier decomposition such that perturbation quanti- 
ties are the real part of a sum of terms oc eico{imip~ imQ.pt), 
m = 0, 1, 2...., one obtains th e following system of equations 
jGoldreich fc Tremainel[l979r ) for the corresponding Fourier 
coefficients, 



iav' — 2Q.u' 



dw' m + *i4 + fmc = 



dr 



dr 



iau m + 2Bv m H 1 = 0, 



ijrWtn , dv'm I (i a ) v ' m I 
c2 dr r 



r 

imuL 



= 0, 



(12) 
(13) 
(14) 



where primes denote perturbation quantities, and the sub- 
script m, being the azimuthal mode number, indicates the 
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Figure 2. Real (solid) and imaginary (dotted) part of the en- 
thalpy perturbation W , for an isothermal, constant surface den- 
sity disc with hp = 10 -3 / 2 and a softening parameter as used in 
KP93 (b = 10 — 4 ) plotted as a function of radius in units of r p . 



m-th Fourier coefficient. Here the velocity perturbation is 
v m = ( v 'm,,u'm), W m = c s Yi' m /Yi is the enthalpy perturba- 
tion, a = m(Q — fi P ), B is the second Oort constant, and <&' m 
is the m-th Fourier component of $ p , being a real quantity. 
The term proportional to /3 is not present in the equations 
of KP93, because they considered a strictly barotropic equa- 
tion of state. Eliminating u m , we obtain a system of ordinary 
differential equations (KP93): 



dr 



171 " \ rrrl 



v (TV 

dWL _ ■ 
dr 

2mn(W m + <S>' m ) 



2mB \ v 

r 



+ 



a 



or 



dr 



(15) 



(16) 



where k 2 = 4B57 is the square of the epicyclic frequency. 

We have solved equations (|15[) and (|16[) using a sixth 
order Runge-Kutta method and outgoing wave boundary 
conditions (see KP93) . In Fig. we show the resulting W m 
for m = 10, for an isothermal, constant surface density disc 
with hp = 10 _3//2 and a very low value of the softening 
parameter b = 10 -4 . The same case was shown in KP93 
(their figure 2), and the results agree very well. 

The imaginary part of W' m is directly related to the 
torque density: 



dT m , ^ lm(WL 

— — = 7rm<I> m E — 

dr ci 



(17) 



and the total torque on the planet r m can be found by in - 
tegrating over the whole disc. iGoldreich fc Tremaind (|l979T l 
showed that the Lindblad torque is carried away by density 
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Figure 3. Total torque on the planet, obtained from linear cal- 
culations, as a function of the softening parameter b. The sur- 
face density is proportional to r~ 3 / 2 (a = 3/2), which means 
that the corotation torque vanishes in the barotropic case. The 
disc is either fully isothermal with h p = 0.05 (dashed line) or 
has a constant aspect ratio h = 0.05 (solid line). Calculations 
similar to those of KP93 with constant h are indicated by the 
dash-dotted line. Also shown are the results of linear calculations 
in 2D by KP93 (for constant h) and the isothermal results of 
TTW02, for 2D as well as 3D. Both 2D studies used a very small 
(KP93) or vanishing (TTW02) softening length. Filled circles de- 
note the torques measured from hydrodynamical simulations for 
1.26 ■ 10~ 5 planet in a disc with constant h, and diamonds 
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the results for isothermal simulations with h v 
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waves, resulting in an angular momentum flux 



p - 



mirrY, 



{Ke{W m ) + C) 



^ m 

dXm{W' m 



dr 



(18) 



Therefore, the Lindblad torque is given by AF m — 
-FmfVout) — Fm(r-ux), where n n and r ou t denote the inner and 
the outer radius of the disc, respectively. We can therefore 
calculate the corotation torque as (KP93): 



AF„ 



(19) 



and the total torques can be found by summation over all 
m . 

4.2 Results 

We start by considering a disc with constant specific vortic- 
ity, which means that the corotation torque vanishes. This 
allows us to look in some more detail at the Lindblad torque 
alone. In Fig. [3] we show the total torqu^B for three differ- 
ent cases: strictly isothermal (dashed line), locally isother- 

1 In all figures, the torque is given in units of S p rpf2p, and is 
divided by q 2 to make it independent of the mass of the planet. 



Figure 4. Lindblad (dash-dotted curves) and corotation (dashed 
curves) torques obtained from an isothermal (h p = 0.05) linear 
calculation as a function of softening parameter b for two different 
density profiles. Black curves: a = 1/2, grey curves: a = —1/2. 
The solid curves denote the total torques, and the horizontal dot- 
ted lines the 3D results of TTW02. 



mal (h constant; solid line) and locally isothermal without 
the term proportional to f3 in equation (|12p (dash-dotted 
line). The latter case is similar to the one considered in 
KP93, and approaches the result of KP93 for small soft- 
ening. Similarly, the strictly isothermal case approaches the 
2D result of TTW02 for small softening. For appropriate 
values of b ~ h, differences between the strictly and locally 
isothermal models can be as large as 30%. In the remain- 
der of this paper, we will restrict ourselves to the strictly 
isothermal equation of state, since the h orseshoe dr ag has 
been analysed for barotropic fluids only (Ward 19910. For 
this smoothing length of b = 0.025, corresponding 

to 0.5 h, results in the total torque being equal to the 3D 
calculations of TTW02. 

We now introduce a corotation torque in the calcula- 
tions by considering different values of a. In Fig- SI we show 
the linear Lindblad and corotation torques for a — 1/2 and 
a = —1/2. Note that in the latter case, the surface den- 
sity gradient is positive, which may be unreali stic except for 
speci al locations in the disc (see for example iMasset et al.l 
boOfj ). but it serves as a good example of a case with strong 
corotation torques. For a = 1/2, the corotation torque is 
much smaller than the Lindblad torque, and is increased 
by a factor of 3 going from b = h to 6 = 0, similar to the 
Lindblad torque. Therefore, the corotation torque is a small 
fraction of the Lindblad torque for all values of b considered 
here. That is not the case for a = —1/2, where for small soft- 



2 In Paardekooper & Papaloizou (2008), horseshoe drag was 
studied in adiabatic flows, where entropy (but not vortensity) 
is conserved. For a locally isothermal equation of state, neither of 
these is conserved, with the consequence that it is unclear what 
the horseshoe drag is in this case. 
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Figure 5. Corotation torque, obtained from a linear calculation 
with a = 1/2 and h p = 0.025 (solid line), h p = 0.05 (dashed 
line), and h p = 0.1 (dash-dotted line). The dotted line gives the 
prediction of equation 1241 valid for h <C b <C 1. 



ening parameters the corotation torque is almost the same 
as the Lindblad torque. The sign of the total torque will 
change for a < —1/2, which is expected from the 2D results 
of KP93 and TTW02. For 3D calculations, the situation is 
different (TTW02) , which is illustrated by the results in Fig. 
[4] for larger softening. For all values of a considered in Figs. 
[3] and 0] we find that for a softening length of 6 — 0.5 h we 
can match our 2D calculations with the 3D work of TTW02. 




0.05 



Figure 6. Corotation torque, obtained from a linear calculation 
with a = 1/2 and h p = 0.05 (solid line), together with the horse- 
shoe orbit drag (dashed line) for our simple estimate of the width 
of the horseshoe region. Also shown is the 3D result of TTW02 
(dash-dotted line). 



from large m ~ 1/6, we replace the sum by an integral, thus 



mKo(mb) 2 



xK (x) 2 dx = \ir 



(23) 



4.3 The limit of large softening 

We now relate our numerical results to the predictions made 
using the torque formula (O given bv lGoldreich fc Tremainel 
(1 19791 ). This is applicable to the situation where the soft- 
ening parameter is significantly larger than the scale height 
and we use it to evaluate the corotation torque in that limit. 
To do this we adopt the a pprox imation for <&' m introduced 
bv lGoldreich fc Tremainel l|l980l ) for m > in the form 



Thus we obtain 



4>p cos{mtp)dip 



2 GM C 



K (mC), 



where Kq is the standard Bessel function and 



(r - r p ) 2 

r 2 
' p 



+ b 2 . 



(20) 



(21) 



This is valid in the important domain of interest where m ~ 
1/6, and r p — t- <~ br p . Summing the torques given by ^ 
over m, using 1|20[1 we obtain 



(22) 



where E p denotes the unperturbed surface density at the 
planets location. Noting that 6 is small and that the domi- 
nant contribution to the sum on the right hand side comes 



E 



r c ,i 



r 4 fi 2 



(24) 



We remark that IWardl (|1992| ) obtained a corresponding ex- 
pression expression for a vertically averaged potential, that 
has the same scaling when 6 is replaced by h p in the above 
equation. 

A comparison of the torques given by equations ([5]) and 
(|24|) indicates that the effect of softening should be signifi- 
cant for 6 >~ lAh p . But note that in addition to requiring 
6/ft > 1 (|24p also formally requires 6 -C 1. This means, as 
we shall see below, that the limit where (|24p applies requires 
rather small h <~ 0.025. 

We now directly compare linear calculations with the 
prediction of equation (I24|l . The results are illustrated in 
Fig. [5] for three different values of h p . Note that in the 
derivation of equation (|24[1 it was assumed that 6 <C 1, so 
we expect the torque to be given by equation (|24[) for h <C 
6<C 1. From Fig. [5] we see that for h p = 0.025 we can nicely 
reproduce equation (|24|) for 0.05 < 6 < 0.2. However, as 
/ip increases there are increasing deviations from equation 
((24)) . For h p = 0.05 and 0.1 < 6 < 0.2, the maximum amd 
minimum deviations are by a factor of 2 and 1.5 respectively, 
while for h p — 0.1 the deviation always exceeds a factor of 
two. 
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5 A COMPARISON OF THE EXPECTED 
HORSESHOE DRAG WITH THE LINEAR 
COROTATION TORQUE 

In this section we compare corotation torques obtained from 
linear calculations to the expected horseshoe drag. We ob- 
tain the latter from equation <(9j which requires an estimate 
of x s which we obtain from equation (110p which simulations 
have indicated gives a good estimate for b ~ h p . We shall 
obtain estimates for the horseshoe drag directly from simu- 
lations and make further comparisons below. 

In Fig. [SJ we show the linear corotation torque (solid 
line), together with the horseshoe drag, obtained as indi- 
cated above, (dashed line) for a disc with a — 1/2, for differ- 
ent values of the smoothing parameter b. For reasonable val- 
ues of 6, the horseshoe drag is stronger than the linear coro- 
tation torque making the torque on the disc more negative. 
Thus the horseshoe drag on the planet is positive and larger 
than the linear corotation torque. For smoothing parame- 
ters 0.02 < b < 0.03, this non-linear torque is a factor 2 — 3 
larger than the linear torque. For smaller values of b, equa- 
tion (I10|) predicts a value of x s that is too large and accord- 
ingly a horseshoe drag that is too large. In reality we expect 
both the horseshoe width and therefore the horseshoe drag 
to reach limiting values for small b and these values should 
be larger than those we obtain here using values of b for 
which equation (|10p app lies. A more complete discu s sion o f 
these aspects is given in lPaardekooper fc Papaloizou! (2009). 
Also shown in Fig. [5] is the corotation torque obtained for 
the corresponding three-dimensional disc (TTW02) which 
has zero softening. For b = 0.02, we can match our 2D lin- 
ear calculation to the 3D result of TTW02. However, for the 
same softening the horseshoe drag is three times as large. 
For a large softening with b — 0.7h p — 0.035, the horseshoe 
drag is equal to the 3D unsoftened linear corotation torque. 
However, numerical results show that for this smoothing 
length, equation (|10[) actually slightly underestimates the 
true value of x s , with the consequence that the intersection 
of the horseshoe drag with the 3D linear corotation torque 
occurs for approximately b = 0.8h p . Lacking a full 3D model 
of the horseshoe region, it is difficult to say what value of 
b to choose so that 2D results match 3D results. However, 
numerical results su ggest that b < 0.02 may be appropriate 
l|Masset et al.|[200r3 ). We come back to this issue in Section 
[7] but it is good to keep in mind that the effects discussed 
in the next sections may well be stronger in 3D calculations. 



6 HYDRODYNAMICAL SIMULATIONS 

In the previous section, we have established that horseshoe 
drag is potentially much stronger than the linear corotation 
torque. We now turn to numerical hydrodynamic simula- 
tions to show that indeed strong deviations from linear the- 
ory are encountered in practice. W e use the RODEO method 
(Paardekooper fc Mellem a 2006b) in two spatial dimensions, 
on a regular grid extending from r — 0.5r p to r = 1.8r p and 
which covers the whole 2-7T in azimuth. Since we want to re- 
solve the horseshoe region for even the smallest planets we 
consider, a relatively high resolution is required: 1024 cells in 
the radial and 4096 cells in the azimuthal direction, making 
the resolution at the location of the planet approximately 
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Figure 7. Total torque on the planet for an isothermal, inviscid 
h p = 0.05 disc with b = 0.03 and q = 1.26 ■ 10" 5 , for different sur- 
face density profiles. The dotted lines indicate the linear torque, 
for increasing a from top to bottom. 

0.0015r p in both directions. Tests have shown that this res- 
olution is sufficient to capture the horseshoe dynamics for 
x s > 0.004. We include all disc material in the torque calcu- 
lation (not excluding any material close to the planet), and 
the planet is introduced with its full mass at t — 0. For low 
mass planets, this does not affect the results. 

First, in Section 16.11 we study the development of the 
horseshoe drag and its relationship to linear theory. Then, in 
Section [6.21 we study the long-term behaviour (saturation) 
of the corotation torque. 

6.1 Development of the non-linear torque 

We start by considering a planet with a relatively high mass, 
q — 1.26 ■ 10~ 5 , which corresponds to 4 Mq around a So- 
lar mass star. Although this is the planet with the largest 
mass we consi der, it has been exp ected to be well within the 
linear regime IjMasset et al. 2006). In Fig. H 

we show that 

for a disc with constant specific vorticity, we can reproduce 
the expected linear torque for all reasonable values of b. For 
very small values of b, one may expect a departure from lin- 
ear theory, since at this point, an envelope may form that is 
gravitationally bound to the planet, which probably should 
be excluded from the torque calculation. We do not con- 
sider this regime here, since we expect to adopt a value of b 
comparable to h in order to account for vertical averaging. 
The important point is that we can match our linear calcu- 
lations to the hydrodynamical simulations in the absence of 
corotation torques (linear as well as non-linear). 

This is further illustrated in Fig. [7] where we show the 
time evolution of the total torque on the planet for different 
surface density profiles. We see that the case with a = 3/2 
nicely falls on the corresponding linear result. Note also that 
this torque is set-up in approximately one orbital period. 
This is to be expected for both the Lindblad and the linear 
corotation torque (see the Appendix). 
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Figure 8. Total torque on the planet in an isothermal h p = 0.05 
disc with a = -1/2 for b = 0.03 and g = 1.26 ■ 10~ 5 for dif- 
ferent magnitudes of the viscosity. The horizontal lines indicate, 
from bottom to top, the total linear torque, the Lindblad torque 
with the horseshoe torque (found using the procedure outlined in 
Section [5J added, and the total linear torque with the horseshoe 
torque added. 




Figure 9. Total torque on a q = 1.26 ■ 10~ 5 planet embedded 
in an isothermal h p = 0.05 inviscid disc for b = 0.03 after 30 
orbits. The solid line indicates the linear torque, and the dashed 
line shows the total linear torque plus the estimated horseshoe 
drag. The dotted line indicates the linear Lindblad torque plus 
the non-linear horseshoe drag. Symbols denote results obtained 
from hydrodynamical simulations for different values of the sur- 
face density power law index a. 



Different surface density profiles give remarkably differ- 
ent results. All cases except a = 3/2 show a departure from 
linear theory, the sign of which is dictated by the vorten- 
sity gradient. This indicates that the corotation torque is 
enhanced with respect to its linear value. This enhancement 
takes approximately 20 orbits to develop after which the 
torques attain steady values for the remainder of the sim- 
ulations. This time period, as we show below, can be un- 
derstood as being a fraction of the libration time scale, and 
is therefore related to a non-linear effect. Note that in all 
cases, linear theory is only valid at early times (less than 
about two orbits) during which as expected the linear coro- 
tation torque is set up on a dynamical time scale. At later 
times, it gets replaced by the non-linear horseshoe drag. 

We take a closer look at the a = —1/2 case illustrated 
in Fig. [5] The horizontal dotted lines indicate, from bottom 
to top, the linear torque, the linear Lindblad torque with the 
horseshoe drag added, and the total linear torque with the 
horseshoe drag added. From the solid line, we see that the 
total torque can be understood as the linear torque, with 
the linear corotation torque replaced by the horseshoe drag. 
We have found this to be true for all values of a considered 
(see Fig. [9j . It is also expected from the discussion given in 
Section [3] 

The magnitude of the non-linear torque depends on the 
magnitude of the viscosity in the disc. Recall that when- 
ever we include viscosity in the model, we do not allow for 
large scale mass flow with respect to the planet (see also 
Section [7|). It is easy to understand this dependence, since 
the horseshoe drag model hinges on vortensity conservation 
while the fluid executes a horseshoe turn. For strong enough 
viscosities, we can recover the linear torque, but note that 
the large values v >~ 10 -3 required correspond to a very 



large a viscosity parameter of a v isc = 0.4. This is because 
only then can the viscosity directly affect the horseshoe turn. 
In Section f6. 21 we will see that lower viscosities are in fact 
sufficient to keep the torque unsaturated. 

Returning to inviscid discs, we show in Fig. [9] the to- 
tal torque on the disc for different values of a, confirming 
that the horseshoe drag indeed replaces the linear corota- 
tion torque. Note that we expect a torque reversal around 
a = — 1, while the linear calculations predict this would only 
happen around a — —3. Note also that steep density profiles 
result in an acceleration of inward migration with respect to 
the linear estimate. 

Linear theory predicts that the torque should scale as 
q 2 . Since in our simple model, the horseshoe width scales as 
g 1 / 2 , the non-linear contribution to the torque (the horse- 
shoe drag) also scales as q 2 . Therefore, it could be misinter- 
preted as a linear effect. However, the time scale on which 
the horseshoe drag term is set up depends on the mass of 
the planet. We will see below that is takes a fixed fraction 
of a libration time scale. This means that although the full 
non- linear torque scales as q 2 , this will not be the case at 
intermediate stages. This is illustrated in Fig. 1101 where we 
show the time evolution of the torque for different mass ra- 
tios. For linear torques, all curves would fall on top of each 
other. This is indeed true at early times, when the horse- 
shoe drag has not yet developed. When the horseshoe drag 
takes over, however, planets of different mass give different 
results; lower-mass planets take longer to develop the non- 
linear torque. 

If we rescale the time axis to the libration time (see Fig. 
ITT)) . the curves fall on top of each other again. This is because 
the time scale to set up the horseshoe drag is a fixed fraction 
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Figure 10. Total torque on a planet in an isothermal h p = 0.05 
disc for b = 0.03 and a = —1/2, for different values of q. The 
dotted line indicates the linear torque. 
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Figure 12. Total torque on a planet in an isothermal h p = 0.05 
disc for q = 1.26 ■ 10 -5 , b = 0.02 and a = -1/2, for different 
values of the viscosity v. 
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Figure 11. Total torque on the planet embedded in an isothermal 
hp = 0.05 disc for b = 0.03 and a = —1/2, for different values of 
q. The time is in units of Tm,, which scales as q~ 1 ' 2 . 
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Figure 13. Total torque on a planet in an isothermal h p = 0.05 
disc for q = 1.26 TO - 5 and a = — 1, for different values of softening 
parameter b and viscosity v. 



executed a horseshoe turn at constant entropy (replacing 
specific vorticity applicable to this case). 



of the libration time (approximately 15%, according to Fig. 
Ill[) . One might expect that, as the scale of the region con- 
tributing the torque is of order H, the time required would 
be on the order of a factor h p smaller than the libration time. 
Thus the results presented here are consistent with this time 
being ab out 2h p Tub- This type of phenomeno n was also dis- 
cussed in iPaardekooper fc Papaloizoul (|200Sl ). where it was 
shown that the development of the non linear torque is due 
to an high density ridge resulting from material that has 



6.2 Long-term evolution 

Lindblad torques give rise to waves that carry away angular 
momentum. Therefore, the planet can continue to exchange 
angular momentum with the disc at Lindblad resonances. 
However, there is no wave transport in the corotation region, 
and therefore in the absence of any other form of transport, 
only a finite amount of angular momentum exchange with 
the planet can occur before the structure of the coorbital 
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region is significantly modified. In other words, the corota- 
tion region is a closed system unless there is some diffusive 
process operating in the disc that can transport angular mo- 
mentum there. In the abs ence of such diffusion, t he corota- 
tion torque will saturate (|Ogilvie fc Lubow| [2003'). 

Saturation is an inherently non-linear process, as it de- 
pends again on a finite width of the corotation region. The 
time scale on which saturation operates is the libration time 
scale, and diffusion must operate on a smaller time scale in 
order to prevent saturation. For a viscosity coefficient v, one 
needs 

v > fr 2 p £l p , (25) 

fsee lMassetJl2002l , who considered this process for a disc with 
constant surface density) . For a level of viscosity that has be- 
come standard in disc-planet interaction studies, v = 10 -5 
all protoplanets considered here are of small enough 
mass that we expect the corotation torque to be unsatu- 
rated. We consider a planet of q — 1.26 -10 , and use a soft- 
ening parameter b = 0.02 to make the libration time scale as 
short as possible to ease the computational burden. In Fig. 
1121 we show the long-term evolution of the torque for three 
different levels of viscosity. The inviscid case shows strong 
libration cycles, before the torque starts to settle to a value 
that is close to the linear Lindblad torque (Fl ~ — lOOOg 2 , 
see Fig. 0}. 

A small viscosity of v = 10~ 7 gives less prominent 
libration cycles, and the torque settles at a value of ap- 
proximately —850. In this case, corotation torques are par- 
tially saturate d. This is in quantitative agreement with 
the analysis of Massct (2001), where it is argued that the 
horseshoe drag should be multiplied by a factor T{z) with 
z — x s (2tyu)~ 1 ^ 3 to account for saturatio n. Choosing the 
smallest o ption for T giv en in lMassell l|200ll) . which gave the 
best fit in lMassetl l|2002l '). we find T = 0.114. This predicts 
a corotation torque of To = 125, while the difference be- 
tween the Lindblad torque and the total torque as measured 
from the simulations (the asymptotic difference between the 
dashed and solid curve in Fig. [T2)l is 150. 

For v — 10~ 5 the libration cycles disappear, indicat- 
ing that the torque remains unsaturated. There is basically 
no evolution of the torque once the non linear contribu- 
tion ha s been set up. This again agrees with the analysis of 
iMassetl (|200lh . where for this value of v we expect T close 
to unity. Note, however, that the maximum torque that can 
be reached is reduced for this relatively high viscosity (see 
also Fig.©. 

Finally, we show in Fig. [T3]that, for a softening param- 
eter b — 0.02, low mass planets will feel a sustained positive 
torque when a = — 1, which corresponds to a mildly posi- 
tive surface density gradient. Although a viscosity v — 10~ 5 
reduces the non- linear torque, making the torque negative 
for b = 0.03, it is necessary for the torque to be sustained. 
For a slightly smaller softening parameter, b — 0.02, we find 
sustained outward migration. Note that from Figs. [3] and 
[4] we expect a smoothing around b = 0.025 to reproduce 
3D effects. Full 3D simulations, together with a 3D under- 
standing of horseshoe dynamics, are needed to see how this 
non- linear torque behaves in a three-dimensional setting (see 
also Section [7}. 
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Figure 14. Schematic overview of the different migration 
regimes. Depending on the magnitude of the viscosity, low-mass 
planets will either experience a linear, unsaturated torque, a non- 
linear, unsaturated torque, or a non-linear, saturated torque. 

7 DISCUSSION 

We have shown that non-linear effects commonly occur in 
the coorbital region, even for low-mass planets. In inviscid 
discs, the corotation torque eventually always becomes non- 
linear, saturating after a few libration cycles, so that only 
the (linear) Lindblad torque survives. Viscosity can prevent 
the torque from saturating, and for a large enough viscosity, 
which however decreases with protoplanet mass, the linear 
torque will be restored. In Fig. Q3] we present a schematic 
overview of the three possibilities in the (log q, log v) plane. 
Two lines define the boundaries between the three regions. 
The boundary between the region where the non linear 
torque is maintained by viscosity and the saturated regime is 
given approximately by the condition that viscous diffusion 
across the coorbital region occurs in one libration cycle. The 
line separating the sustained horseshoe drag regime from the 
linear regime is expected to occur where viscous effects be- 
come large enough to disrupt horseshoe turns. 

Note that the regime for which the corotation torque is 
linear, occupies only a small fraction of the parameter space. 
Note also that the borders between the different regions are 
not razor-sharp. In reality the non-linear torque can be par- 
tially saturated, or be reduced by viscosity. However, the 
overall picture is clear: the corotation torque is almost al- 
ways non-linear. 

For density profiles with close to a = 3/2, the non- 
linear corotation torque, or horseshoe drag, does not play 
a major role. In an isothermal disc with a = 1/2, a com- 
monly adopted value for a locally isothermal disc, the devi- 
ation from the linear torque can be up to 30% in the invis- 
cid case (see Fig. [JJ. If a relatively large viscosity is used, 
such non-linear be haviour can be markedly reduced ( see also 
the discussion in iPaardekooper fc Papaloizoul 12008! ) . Note 
also that for a given viscosity, linearity is always restored 
for small enough protoplanet masses (see Fig. I14[) . Large 
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viscosities, numerical and/or imposed, in combination with 
inadequate numerical resolution for representing horseshoe 
turns for low-mass planets, conspire to make the non-linear 
behaviour described here less apparent. 

iMasset et~ al. (2006) reported non- linear behaviour for 
planets of higher mass than considered here, also claimed 
to be due to the action of the horseshoe drag. However, 
since the horseshoe width x B was enhanced relative to its 
value for low-mass planets, th e non-linear torque was much 
stronger. IMasset et all <|2006f > define non-linearity as a de- 
parture from the torque being proportional to q 2 . In the light 
of our findings, however, it should be noted that all planets 
in fact show non- linear behaviour in the inviscid limit. For 
low-mass planets, this non-linearity is less obvious because 
the horseshoe drag is also proportional to q 2 . 

For moderately positive density gradients (a —1), the 
non-linear torque is strong enough to reverse the sign of the 
total torque. In realistic discs, positive gradients will prob- 
ably only be realised at special locat i ons, b ut note that, 
unlike as was argued in IMasset et al. I (|2006h . the gradient 
does not have to be extremely sharp to stop inward migra- 
tion. Almost any positive surface density gradient can act 
as a 'protoplanet trap'. 

Barotropic discs with a < can also serve as a model 
for what happens when strong corotation torques arise that 
are not due to the radial vortensity gradient. In adiabatic 
discs, for example, there is a strong contribution from any 
radial entropy gradient (IPaardekooper fe Mellemal 2008 ; 



weighted average (see lMassetl \200i )): 

4 S™ 00 P{rp,<f P ,z)xj{b = z)dz 



(26) 



IPaardekooper fc Papaloizou 20081 ) . iBaruteau fc Massetl 



However, this procedure requires a proper estimate for x s 
for 6^0, where equation (| 10p cannot be used. This issue is 
discussed in a paper. 



8 SUMMARY AND CONCLUSION 

We have analysed the corotation torque on an embedded 
planet in a barotropic disc, and shown that this torque is 
non-linear in general after a few orbits. The linear corota- 
tion torque, which is set up on a dynamical time scale, is 
replaced by horseshoe drag, which is stronger in all cases we 
have considered. This process completes in approximately 
one tenth of the libration time scale nib or ~ 2/i p n;b in our 
case. 

For discs with large vortensity gradients, a strong de- 
parture from linear theory is observed. We have shown that 
in discs with moderately positive density gradients (E ex r) 
non-linear effects can reverse the total torque, leading to 
outward migration. In particular, this may occur without a 
very abrupt strong surface density transition of the typ e con- 
sidered by ('protoplanetary trap'. M asset et alj r2006) . and 
may lead to the halting of the inward migration of low-mass 
planets. 

After one libration time, saturation sets in unless some 
form of viscosity is able to restore the original density pro- 



(200$) showed that linear theory predicts a contribution file within nib- For the standard value of v 



10" 



from a radial entropy gradient, and argued that this would 
be strong enough to rever s e the total torque. However, 
IPaardekooper fc Papaloizou! (|2008h showed that this linear 
contribution is in fact small, and that it is a non-linear 
effect that changes the s ig n of t he torque as observed in 
Paardekooper fc Mellemal (|200 6a) . and more recently in 
Kiev fc Cridal (|2008l ). The non-linear contribution arises in 



a similar manner to that discuss ed in this paper (see also 
IPaardekooper fc Papaloizoull2008l ). 

We have greatly simplified the problem by keeping the 
planet on a fixed orbit and choosing the viscosity law such 
that there is no accretion flow. This way, there is no radial 
mass flow with respect to the planet. Material that flows 
past the planet from the inner to the outer disc (or the other 
way round) exerts an additional torque on the planet (see 
M asset 2001). The effect of such a radial flow is to introduce 
an asymmetry in the horse shoe region, which lies at the 
basis of Type III migration j|Masset fc Papaloizou! |2003 ). It 
remains to be investigated how these processes affect the 
current analysis. 

IMasset et ail (|2006l ) noted that the effects of non- 
linearity were much stronger in 3D, probably because the 
width of the horseshoe region is larger compared to 2D sim- 
ulations that necessarily have a softening parameter of the 
order of the scale height of the disc. The full 3D structure of 
the horseshoe region remains to be investigated, but if the 
velocity field is essentially two-dimensional, one could regard 
the vertical structure of the horseshoe region as stacked lay- 
ers of 2D horseshoes, with a smo othing equal to the vertical 
distance to the midplane (see also lMasset et al" ] |2006h . Then, 
the term xf in equation (|9]) would be replaced by a density- 



sponding to a v isc = 0.004, the corotation torque on low-mass 
planets in the Earth mass range can be sustained by such 
action. 
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APPENDIX 

TIME DEVELOPMENT OF THE LINEAR 
COROTATION TORQUE 

We here consider the time development of the linear corota- 
tion torque consequent on the insertion of a perturbing pro- 
toplanet potential into an unperturbed disc. We specialise 
to the limit of a low mass planet and adopt a local Cartesian 
coordinate system (x,y) with origin at the centre of mass of 
the planet, and the x axis pointing radially outwards. The 
system rotates uniformly with the Keplerian angular veloc- 
ity, flp, at the orbital location of the planet. In this frame 
the differential rotation of the disc is manifest through a 
linear shear v = (v x ,v y ) = (0, — 3£l p x/2). Thus the back- 
ground vorticity is constant and the vortensity gradient is 
proportional to the gradient of 1/E. 

The governing equations are the inviscid form of the 
basic equations (fTJ and . As it builds up from nothing, the 
disc response has to be linear for some period of time after 
insertion of the protoplanet so we linearise the governing 
equations treating $ p as a linear perturbation. The linear 
response approaches its final value in a characteristic time 
independent of the protoplanet mass, in fact we shall see that 
this is characteristically the orbital time scale. Accordingly 
a full linear response is guaranteed for a sufficiently small 
protoplanet mass. 

For simplicity we initially discuss the large softening 
case for which the pressure response is small (see Section 
14.3) 1 . We then go on to discuss the more general case. 



Calculation of the linear response 

The linearised forms of the components of ((2)1 in the local 
system are 



Dv' x , <9$Gp 

__ 2 ^ = ___, 



and 



Dv 'v ! n ' 

— — — H — S2n« x 

Dt 2 P 



gggp 
dy 



(27) 



(28) 



where perturbation quantities are denoted with a prime and 
the convective derivative applies to the unperturbed flow 
such that 



D_ 

Dt 



d 3Q p x d 
di 



2 dy' 



(29) 



These are the local analogues of the global linear equations 
(|12[) and (|13|) . $Gp here denotes the generalised potential 
which may be taken to be the sum of the perturbing pro- 
toplanet potential and the enthalpy perturbation. For the 
large softening model, $Gp may be taken to be the proto- 
planet perturbing potential alone. Note we have not made a 
Fourier decomposition at this point. 

We find it convenient to work with the Lagrangian 
displacement £ = (£x>£y), which is such that (v' x ,v' y ) = 
(D£ x / Dt, D£ y I Dt + 3Q p £ x /2). In terms of this (J27J) and 
take the form 



D 2 i 



on £5l _ qoV 
Dt P Dt p ^ 



&tv , 90 D & _„3 $ Gp 

Dt 2 + p Dt ~ dy ' 



Gp 



dx 



and 



(30) 



(31) 



We now perform a Fourier decomposition of $g p in y as- 
suming periodicity on a length scale L y . 



®G P = Tie 2j bn ( x ' 1 ) ex P( ink yoy) , 



(32) 



where k y o = 2n/L y and Tie denotes that the real part is 
to be taken. We remark that to relate to a full cylindrical 
annulus, L y — > 27rr p , y — > r p tp, and n — * m, where m is the 
usual azimuthal mode number. Using this we obtain 



2%. 2C1 ^--3Q 2 f 
Dt p Dt p4 * 



^ ^7 exp(ink y0 y), and 
(33) 



D 2 £ D£ x °° 

v + 2D.p——^- = — ^2 ink y0 b n exp(ink y0 y). 



Dt 2 



(34) 



In order to solve these equations we make the approxi- 
mation of neglecting D 2 £ x /Dt 2 . This is a common approxi- 
mation that is made when considering corotation/horseshoe 
dynamics. It has the effect of removing epicyclic oscillations, 
and therefore Lindblad torques, but as we shall see it does 
not interfere with the corotation torque. We then use ([33} 
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Figure 15. Linear density response for a 



-2, b = 0.03, h v 



0.05 and q = 1.26 ■ 10 — J . The color scale has been adjusted to 
highlight the density ridge at corotation, which gives rise to the 
linear corotation torque. The planet is located at r = 1 and <f> = it. 



to eliminate £ x from (|34p and thus obtain a second order 
equation for £ y : 



1D% 
3 Dt 2 

where 



= - c» exp(ink y0 y), 

n=0 



2 a 2 fe n 
3fi p dxdt 



ink y o(b„ — xdbn/dx). 



(35) 



(36) 



This may be integrated to obtain £ y and then £ x found from 
the used approximate form of (|33f > - The integration process 
is aided by noting that y + 3Q, p xt/2 is invariant under D. 
Note too that in order to ensure that the boundary condi- 
tion, that all disturbances vanish at t = 0, is satisfied, an 
appropriate function of y + 3Q p xt/2, which plays the role 
of an integration constant, may be added (D/Dt operating 
on such a function will be zero). Following the above proce- 
dures, which guarantee there will be no singularities in the 
solution, we obtain 



= -3 exp(ink y0 (y + 3n p xt/2)) x 



c n exp(— SinkyoQpXt' /2))dt' , (37) 



= - + m £ iZ eMmkyoy) 

p n = 



(38) 



/' 

Jo 



, exp(— SinkyoClpSct /2))(t — t')dt' . 



(39) 



We note that in the above and other similar integrands, 
unless otherwise indicated, quantities are evaluated at the 
time t' . We may now find the density perturbation (see Fig. 
[15] for and example) from 

E' 



-V • ££ 



(40) 



and then evaluate the corotation torque acting on the planet, 
r c ,Un) by doing the torque integral 



/<9<F 



Gp 



dxdy, 



(41) 



where we have adopted a multiplicative factor equal to an 
orbital radius r p in order to convert a force in the y di- 
rection into a torque. Substituting the Fourier expansion of 
d<&Gp/dy and integrating over y we find, after some elemen- 
tary algebra, that r C] u n = T a + Th, where 



r -^T 



^^^-n 2 b n cos Xn(x,t')dt' dx, (42) 



with Xn{x,t') — 3nko y tt p x(t — t')/2, and 



27rr p \ - 



dx 



-nb n sin Xn(x,t')dt'dx, (43) 



with d n — b n — xdbn/dx, and /„ = (2/3Q p )d 2 b n / dxdt' . 

The expression for Fb can be integrated by parts with 
respect to t' and combined with that for T a to yield the 
following expression for the total torque: 



c,lin 



2tt 



E 



dx 



n 2 b n cos Xn(x, t')dt'dx. 



(44) 



Equation (|44|l shows how the corotation torque devel- 
ops with time after a protoplanet is inserted at t = 0. We 
shall see that this is on the time scale of the orbital period. 
Consider first the case when the b n are independent of time 
as would be the case for a large softening length compared to 
the scale height where there is negligible pressure response. 
Then the integral with respect to t! can be performed with 
the result that 



47rr p f°° ^ 

6U P J-™ n=0 



ng n (x = (L y / (3nnQ p t))^^-d(, 



(45) 

where g n (x) = b„(x) 2 (d'E/dx). 

We remark that most of the contribution to the inte- 
gral comes from f ~ 1. The corresponding value of x ~ 
L y / [3n-nQ.pt) ~ 2r p /(3nf2 p t). Thus at early times the coro- 
tation torque has contributions from mainly large radii but 
as time progresses the contributing region contracts towards 
the corotation circle. To estimate the time involved, we note 
that for large softening, the value of n giving the dominant 



Fourier components for the corotation torque is expected to 
be n ~ 1/6, with contributions from larger values of n being 
reduced because of smoothing. We then see that the region 
contributing to the corotation torque is within 6r p of the 
corotation circle within an orbital period independently of 
its size. Thus the linear corotation torque in this model is 
established in about an orbital period. 

The limiting value is easily obtained by letting t — > 
oo in the above expression. The result after performing the 
integration is 
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r c ,iin = -J^ J2 nb n (0) 2 (dZ/dx) x=0 . (46) 

6U P n=0 

This is e xactly what is obtaine d from the expression 
given by iGoldreich fc Tremainel |l979 : ) recalling that the 
vorticity in the background local model is constant and equal 
to f2 p /2 making the vortensity gradient proportional to the 
gradient of 1/S. To illustrate the density response in the 
neighbourhood of corotation we plot the linear surface den- 
sity response for the case with a = —2, 6 = 0.03, h p = 0.05 
and q = 1.26 ■ 10~ 5 in figure [TU This was calculated fol- 
lowing the procedure described in section [4] In addition to 
the strong wake produced by the Lindblad torques there is 
both an overdensity leading and an underdensity trailing the 
planet. These narrow features are localized on the corotation 
circle as expected from the analysis presented here. These 
produce a positive torque acting on the planet. 

Although (|46|) has been obtained assuming the soften- 
ing length Is large compared to the scale height, we expect 
it to apply more generally. When the softening parameter is 
not large and the generalised potential, <!>Gp is used, the 
added enthalpy perturbation results in the 6„ being un- 
known fun ctions of time. However , the linear response cal- 
culation of IGoldreich fc Tremainel (|l979l ) indicates that the 
scale of the response is \x\ ~ H, with the dominant val- 
ues of n ~ r p /H. Equation (|44|l then indicates only values 
of \t — t'\ <~ r p /(nHQ p ) ~ fip 1 will contribute signifi- 
cantly to the torque. Thus when <&Gp varies slowly com- 
pared to the orbital period, expected as the linear response 
approaches its steady value, b n may be taken to be locally 
constant in time, and then from th e above discussion, the 
corota tion torque will be the local IGoldreich fc Tremainel 
l| 19791 ) value. The situation is also similar at early times. 
For t <~ r p /(nHQ p ) ~ Q~ , the cosine term in equation 
(I44p may be replaced by unity. This then implies that at 
early times 

47r 2 r p ^ 3nHQ p t I dE , 2 \ 

where the angled brackets denote an appropriate integral 
mean. When t approaches £l p 1 , this becomes the same mean 
of torques of the form (|4j|. 

Note that the time to establish the linear response and 
corotation torque, the latter expected to be finite, does not 
depend on the size of the perturbation, or accordingly, q. 
This time should be an intrinsic time which can only be 
a multiple of the orbital period. Accordingly for sufficiently 
small q we can ensure that the system is in the linear regime 
when the torque is set up so that the linear formulation 
should be valid then. 



